#pragma once
#include <complex>
#include <cmath>
#include "algebra.h"
#include "refc_ptr.h"

// M_PI isn't part of the C99 standard and as a (sad) result also isn't part
// of the C++ standard as such. However, most compiler do still define the
// macro in cmath / math.h; this is here for the compilers that don't define
// the macro.
#ifndef M_PI
#define M_PI 3.14159265358979323846264338327f
#endif

namespace krryn{
namespace audio{
	// this wrapper makes it easier to substitute with a faster implementation

	template<typename _Tx>
	const _Tx sin(const _Tx a_Number){
		return std::sin(a_Number);
	}

	// Direction of the FFT transformation.
	// The sign indicates forward or inverse FFT
	// with -1 begin forward and +1 being inverse
	enum fft_direction_e{
		forward = -1,
		inverse = +1
	};

	// The input of this class (both for bit_reverse_copy and for transform)
	// should be a array of complex numners [real, imaginary] such that there
	// are _Size amount of complex numbers stored in the array. The total size
	// of the array is thus 2 * _Size.
	//
	// The output of this class has the positive end of the frequency spectrum
	// in [0, _Size] and the negative part of the frequency spectrum in 
	// [_Size + 3, 2 * _Size]. The DC is at thus _Size + 1. Again, the data i
	// stored as [real, imaginary] in the output array.
	//
	// Notice that if the input of this class consists of real numbers (eg. 
	// complex numbers that look like [real, 0]) then the output is mirrored
	// in the negative frequency spectrum; which can then be ignored.
	template<typename _Type, fft_direction_e _Direction, size_t _Size>
	class fft_imp{
	public:
		typedef std::complex<_Type> complex_t;

		fft_imp() : m_Next() {}

		// We rely on the compler to be smart enough to remove this function for 
		// all recursive template declarations because, otherwise the executable 
		// size will grow unnecessary :-).
		void bit_reverse_copy(_Type *a_Data) const{
			const size_t l_FullSize = _Size * 2;
			for(size_t i = 0, j = 0; i < l_FullSize; i += 2){
				if(j > i){
					std::swap(a_Data[j], a_Data[i]);
					std::swap(a_Data[j + 1], a_Data[i + 1]);
				}

				size_t l_Mask = _Size;
				for(; l_Mask >= 2 && j >= l_Mask; l_Mask /= 2){
					j -= l_Mask;
				}
				j += l_Mask;
			}
		}

		// This routine expects the even indices of a_Data to contain the real
		// part and the odd indices to contain the imaginary part of the complex
		// numbers.
		void transform(_Type *a_Data) const{
			// We rely on template meta programming to do the recusion for us, 
			// so we don't have to do it at runtime. For more information see:
			// http://ubiety.uwaterloo.ca/~tveldhui/papers/Template-Metaprograms/meta-art.html
			//
			// We also rely on the Danielson-Lanczos Lemma that shows that a 
			// fourier transform of length k is the sum of two other fourier
			// transformations. One of 0 < x < k / 2 and one of k / 2 < x < k.
			// We also know that the identity of the transform is when x = 1, so 
			// we can end the recursion there. (Eg. when the input is 1 we don't 
			// have to perform any operations because the result is 1).
			//
			// We use these facts to perform the recursion.
			m_Next.transform(a_Data);
			m_Next.transform(a_Data + _Size);

			// Notice that, since l_Theta is completely known at compile-time
			// we should also be able to calculate the sine of l_Theta at 
			// compiletime, removing the costly sin() function calls. As shown
			// in the aforementioned article.
			const _Type l_Theta = _Direction * static_cast<_Type>(M_PI) / _Size;
			const _Type l_Temp  = math::sin(l_Theta);

			complex_t l_Prev(-2 * l_Temp * l_Temp, math::sin(l_Theta * 2));
			complex_t l_Root(1, 0);

			for(size_t i = 0; i < _Size; i += 2){
				const complex_t l_Temp(
					a_Data[i + _Size + 1] * l_Root.real() 
						- a_Data[i + _Size] * l_Root.imag(),

					a_Data[i + _Size] * l_Root.real() 
						+ a_Data[i + _Size + 1] * l_Root.imag()
				);

				a_Data[i + _Size + 1]	 = a_Data[i + 1] - l_Temp.real();
				a_Data[i + _Size]		 = a_Data[i] - l_Temp.imag();
				a_Data[i + 1]			+= l_Temp.real();
				a_Data[i]				+= l_Temp.imag();

				// update root of unity
				const complex_t l_Add(
					l_Root.real() * l_Prev.real() 
						- l_Root.imag() * l_Prev.imag(),

					l_Root.imag() * l_Prev.real()
						+ l_Root.real() * l_Prev.imag()
				);

				l_Root += l_Add;
			}
		}
	private:

		const fft_imp<_Type, _Direction, _Size / 2> m_Next;
	};

	// template specialization to indicate the 
	// end of the recursion that is used in the
	// actual fft_imp class (see fft_imp.m_Next
	// for more details).
	template<typename _Tx, fft_direction_e _Direction>
	class fft_imp<_Tx, _Direction, 1>{
	public:
		void transform(_Tx *a_Data) const {}
	};

	// Wrapper / convenience class for fft_imp
	//
	// Notice, if you can't guarantee _Size at compile-time this isn't for you
	// however, if you can guesstimate the size of the array, you can zero-fill
	// the rest of the elements. The same is useful if you can't guarantee that
	// _Size will be a power of two.
	template<typename _Tx, size_t _Size>
	class fft{
	public:
		typedef core::refc_ptr<fft<_Tx, _Size> > ptr;

		fft() : m_Forward(), m_Inverse() {}

		size_t size() const{
			return _Size;
		}

		void transform_forward(_Tx *a_Data) const{
			m_Forward.bit_reverse_copy(a_Data);
			m_Forward.transform(a_Data);
		}

		void transform_inverse(_Tx *a_Data) const{
			m_Inverse.bit_reverse_copy(a_Data);
			m_Inverse.transform(a_Data);

			// normalize the data
			const _Tx l_Reciprocal = 1 / static_cast<_Tx>(_Size);

			for(size_t i = 0; i < _Size; i++){
				a_Data[i] *= l_Reciprocal;
			}
		}
	private:

		const fft_imp<_Tx, forward, _Size> m_Forward;
		const fft_imp<_Tx, inverse, _Size> m_Inverse;
	};
}
}
